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Abstract 

In molecular simulations, one of the most difficult points is to track the real dynamics 
of many-body systems from the first principle. The present study shows that step-size 
dependences have an unexpected effect on simulation results, even when we use the 
standard high-precision numerical integrators to apply to a simple system with a two- 
body interaction. The validity of our analysis is checked by the theory of adiabatic 
approximations . 

Recently, various numerical integration techniques are developed to simulate 
the dynamics of complex systems with many-body interactions. Above all, sym- 
plectic integrators are important because they exactly preserve the symplectic 
form = dp A dq = J2i=i,n dpi A dqi which any natural Hamiltonian H{q, p) 
must have. Thus, it is very important to investigate and estimate how these 
symplectic integrators improve our capacity of tracking many-body dynamics. 
We consider the following simple system to be tested: 

H=lipl+pl + q!ql). (1) 

The equations of motion remain invariant under the transformation (t' — at, q' 
-q,p' — >■ ^p), where t is a time variable and a is an arbitrary real and dimen- 
sionless quantity which does not vanish. Thus, our simulations do not depend on 
the scale of time. It is known that this quartic potential system shows chaotic 
behavior, as is predicted by the non-integrability proof Furthermore, the sys- 
tem (|l]) allows not only chaotic but also intermittent behavior, which significantly 
affects the accuracy of the higher-order symplectic simulations ^. 

In this paper, we investigate the time step-size dependences of its simula- 
tions by two different explicit type higher-order symplectic integrators, namely 
(i) Suzuki's third-order symplectic integrator and (ii) Ruth's third-order sym- 
plectic integrator, which are obtained by the real decompositions of exponential 
operators; 

^{A+B)At ^ ^ciAAt^diBAt^CiAAt^diBAt^csAAt^dsBAt _^ q (^At^) (2) 
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where At is a time step-size and the coefficients Ci, C2, C3, dij c?2, are given as 



ci = d3 = 0.2683300957817599 ■ ■ ■ 
(i^ = C3 = ci + 0.651331427356399 ■ ■ ■ 

C2 = d2 = di- 0.839230460347997 • • • , (3) 

in (i)Suzuki's case[§| or they are given as 

ci = 7/24, di = 2/3, C2 = 3/4, d^ = -2/3, o, = -1/24, 4 = 1 (4) 

in (ii)Ruth's case|^. The main reason why we use third-order symplectic inte- 
grators here hes in the fact that they are already higher-order than the popular 
leap-flog method (the second-order integrator) and it is also guaranteed that there 
are no other third-order symplectic integrators with real coefficients^]. In cases 
of separable Hamiltonian systems H{q,p) = K{p) + V{q), the implementation 
of these third-order symplectic integrators is straightforward as follows: 







-CiAtrg(qr(°)), 




= q^'^ + diAtKpip'-'^] 


p(2) 


= p« 


-C2Atyg(g«), 




= q^'^ + d2AtKpip^^^] 


p(3) 


= p(2) 


-C3Atrg(q(2)), 


g(3) 


= + d-iAtKp{p^^^ 



(5) 

where ■^V{q) = Vq{q), -^K{p) = Kp{p). In our simulations here, we al- 
ways flx the initial condition as (q'l, ^'2,^1,^2) = (1000,0.002,0,0) to cause a 
strong intermittency . Around this initial condition, we can employ the adia- 
batic approximation as 

H=\pI + {\pI + \kql) = E, (6) 

where qf = k ^ Const, is a slow variable. A fast dynamics is simply described 
by the harmonic oscillator Hp = ^pl + |A;g| = Ep with the spring coefficient 
k = qf. Thus, the adiabatic invariant J is given by the formula 

12 12 2 

J = ^El/k = 2P2+^2^iQ2 ^ ^4/10002 = 0.002, (7) 

because k = qf = 1000^ and E'^ = 4 hold at the initial condition. Thus, we have 
a slow dynamics given by the Hamiltonian H = ^pf + J ■ \qi\ and therefore we 
can obtain the analytic solution of this slow dynamics 

q,{t) = q,{t = 0)-^Jt'. (8) 

By considering that the hypothesis of this adiabatic approximation is based on 
the assumption ^ 1, an adiabatic transition is predicted to occur at t = T 
given by 



T = J2qi{t = 0)/J = J2 x 1000/0.002 = 1000, (9) 
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where qi{t = T) ^ 0. Figure 1 shows that the third-order symplectic integra- 
tor with a step-size At = 0.001 can correctly track the adiabatic transition at 
t = 1000. Figure 2 shows that a shght step-size difference cause the great dis- 
crepancy in traclcing the trajectory {qi{t) , q2{t)) , even though each third-order 
symplectic integrator with the step-size around At = 0.001 can give the correct 
behavior of adiabatic transitions like Fig.l. Figure 3.(a)-(c) show that this effect 
of step-size dependences of the dynamical variable q2 at the fixed transition time 
t = 1000 is extraordinary complex beyond our imagination. Remark that this 
effect does not depend on the choice of integrators as is indicated in Fig.3.(a)-(c). 
This result suggests that the following usual argument that "if At = 0.002 is not 
O.K, then At = 0.001 is maybe O.K." is generally untrue unless the step-size At 
is sufficiently small such as < At < 0.0001 like Fig.3.(c). Thus, this kind of 
complex behavior of step-size dependences can be a real challenge to the solid 
progress in the studies of many-body dynamics simulations, since the intractabil- 
ity here is closely related to the non-integrable character of generic many-body 
systems]^, 0, whether the simulations are classical or quantum ones|P]. 
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Figures 



Fig. l:The time evolution of tlie adiabatic invariant J{t) calculated using the 
third-order symplectic integrator for At = 0.001 with the initial condition qi = 
1000, q2 = 0.002 and pi = P2 = 0. An adiabatic transition is observed just at 
t = 1000, as is predicted in Eq. (Q). 
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Fig. 2: The trajectories of {qi,q2) calculated using the third-order symplectic 
integrator for At = O.OOlOl(Green), At = 0.00103 (Blue), At = 0.00105(Red), 
and At = 0.00107(Pink) with the initial condition = 1000, ga = 0.002 and 
Pi = P2 = 0. 
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Fig. 3(a): The dynamical variable q2 at t = 1000 tracked from the unique 
initial condition qi = 1000, q2 = 0.002 and pi = p2 = at t = using two dif- 
ferent third-order symplectic integrators (Suzuki's integrator (Green) and Ruth's 
integrator (Blue)) is plotted against the time step size At (0.001 < A.t < 0.002). 
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Fig. 3(b): The dynamical variable q2 at t = 1000 tracked from the unique 
initial condition qi = 1000, q2 = 0.002 and pi = p2 = a.t t = using two dif- 
ferent third-order symplectic integrators (Suzuki's integrator (Green) and Ruth's 
integrator (Blue)) is plotted against the time step size At (0.0001 < A.t < 0.001). 
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Fig. 3(c): The dynamical variable q2 Bit — 1000 tracked from the unique initial 
condition gi = 1000, 52 = 0.002 and pi — p2 — Bit t — using two different 

third-order symplectic integrators (Suzuki's integrator(Green) and Ruth's inte- 
grator(Bhie)) is plotted against the time step size At (0.00001 < At < 0.0001). 
About iterations are required to estimate each 52 (^ = 1000; At). 
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